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I. INTRODUCTION 



Analytical treatment of problems having quenched disorder is usually difficult. There 
are few models having nontrivial quenched disorder that can be solved exactly. In this 
paper, we obtain exact results for the non-equilibrium properties of the random-field Ising 
model (RFIM) on the Bethe lattice. We consider the single-spin-flip Glauber dynamics of the 
system at zero temperature, as the external magnetic field is slowly varied from — oo to +00. 
As the field increases, the magnetization increases as groups of spins flip up together. This 
model has been proposed as a model of the Barkhausen noise by Sethna et al [T[] (see also 
0). In this paper, we set up the exact self-consistent equations satisfied by the generating 
function of the distribution of avalanche sizes, and analyze these to determine the behavior 
of the avalanche distribution function on the Bethe lattice. 

The study of the equilibrium properties of the RFIM has been an important problem in 
statistical physics for a long time. In 1975, Imry and Ma |J, showed that arbitrarily weak 
disorder destroys long-ranged ferromagnetic order in dimensions d < 2. The persistence of 
ferromagnetism in d — 2 was a matter of a long controversy, but has now been established 
0]. A recent review of earlier work on this model may be found in ||. As far as an exact 
calculation of thermodynamic quantities is concerned, there are only a few results. For 
example, Bruinsma studied the RFIM on a Bethe lattice in the absence of an external field 
and for a bivariate random field distribution [[J. There are no known exact results for the 
average free energy or magnetization, for a continuous distribution of random field, even at 
zero temperature and in zero applied field. 
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FIG. 1. The Hysteresis loop of magnetization M versus the external field h for RFIM. The 
zoomed figure shows the small jumps in magnetization that give rise to the Barkhausen noise 



The non-equilibrium properties of the RFIM has attracted a lot of interest lately, arising 
from the observation by Sethna et al [0 that its zero-temperature dynamics provides a 
simple model for the Barkhausen noise and return point memory. Barkhausen noise is the 
high frequency noise generated due to the small jumps in magnetization observed when 
ferromagnets are placed in oscillating magnetic fields [Fig. ^ . Understanding and reduction 
of this noise is important for the design of many electronic devices J7|. Experimentally it 
is observed P-fl] that the the increase of the magnetization occurs in bursts that span 
over two decades of size and the distribution of burst (avalanche) sizes seems to follow a 
power law in this range. Similar avalanche- like relaxational events are also observed in 



other systems, for example, the stress-induced martensitic growth in some alloys fll2| . This 

as an 



power-law tail in the event-size distribution was interpreted by Cote and Meisel [11 



example of self- organized criticality. But Percovic et al have argued that large bursts are 
exponentially rare, and the approximate power-law tail of the observed distribution comes 



from crossover effects due to nearness of a critical point. Recently Tadic |13| has presented 
some evidence from numerical simulations that the exponents for avalanche distribution can 
vary continuously with disorder. Our results about the behavior of the avalanche distribution 
function also relate to this question whether any fine-tuning of parameters is required to 
see power-law tails in the avalanche distribution in the RFIM, and if the exponents can be 
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varied continuously with disorder. 

The advantage of working on the Bethe lattice is that the usual BBGKY hierarchy of 
equations for correlation functions closes, and one can hope to set up exact self-consistent 
equations for the correlation functions. The fact that Bethe's self-consistent approximation 
becomes exact on the Bethe lattice is useful as it ensures that the approximation will not 
violate any general theorems, e.g. the convexity of thermodynamic functions, sum rules. In 
the presence of disorder, in spite of the closure of the BBGKY hierarchy, the Bethe approxi- 
mation is still very difficult, as the self-consistent equations become functional equations for 
the probability distribution of the effective field. These are not easy to solve, and available 



analytical results in this direction are mostly restricted to one dimension |14| , or to models 
with infinite-ranged interactions [[HJ. On the Bethe lattice, for short-ranged interactions 
with quenched disorder, e.g. in the prototypical case of the ±J random-exchange Ising 
model, the average free energy is trivially determined in the high temperature phase, but 
not in the low-temperature phase. It has not been possible so far to determine even the 
ground-state energy exactly despite several attempts [ITB . 



Calculation of time-dependent or non-equilibrium properties presents its own difficulties, 
even in the absence of disorder. Usually, for d > 1, one has to resort to the limit of coordi- 
nation number becoming large, with interaction strength scaled suitably with coordination 
number to give a nontrivial thermodynamic limit ||17|| . The large-d limit in the self-consistent 
field approximation for quantum- mechanical problems is similar in spirit [|18f| . 

The RFIM model on a Bethe lattice is special in that the zero-temperature nonequilibrium 



response to a slowly varying magnetic field can be determined exactly fll9| . To be precise, 
the average non-equilibrium magnetization in this model can be determined exactly if the 
magnetic field is increased very slowly, from — oo to +00, in the limit of zero temperature. 
It thus provides a good theoretical model to study the slow relaxation to equilibrium in 
glassy systems. The dynamics is governed by the existence of many metastable states, with 
large energy barriers separating different metastable states. We hope that this study of non- 
equilibrium response in this model would help in the more general problem of understanding 



the statistical mechanics of metastable states in glassy systems. 

A brief summary of our results is as follows. We derive the exact self consistent equations 
for the generating function of the avalanche size distribution function Q(x) on the Bethe 
lattice. This is a polynomial equation in Q(x) and x, in which the coefficients depend on 
the external field h, and the distribution of the quenched random fields. We can solve 
these equations explicitly numerically and thus determine the qualitative behavior of the 
distribution of avalanches for any distribution of the quenched random fields. The behavior 
depends on the coordination number z, and on the details of the distribution function. 
We work out the distribution of avalanches explicitly for a rectangular distribution of the 
quenched fields, for the linear chain (z = 2), and the 3-coordinated Bethe lattice. In 
both cases, one finds only exponential decay. We also studied other unimodal continuous 
distributions, e.g. when the random field distribution is gaussian, or of the form Prob(hi) = 
■^sech 2 (^), also for large z. We find that, for z > 4, there is a regime of disorder strengths 
for which the magnetization shows a jump- discontinuity ( "first order transition"), but the 
avalanche distribution, averaged over the hysteresis loop also shows a power law tail of the 
form s~ 5 / 2 ("critical fluctuations"). 

The paper is organized as follows. In section [TJ, we define the model precisely. In 



section pl| , we briefly recapitulate the derivation of self-consistent equations for the magne- 
tization in our model, and then use a similar argument to construct the generating function 
for the avalanche distribution for arbitrary distribution of the quenched random field. We set 
up a self-consistent equation for the probability Q n , that an avalanche propagating in sub- 
tree flips exactly n more spins in the subtree before stopping. The probability distribution 
of avalanches is expressed in terms of this generating function. In section [IV], we consider 
the special case of a rectangular distribution of the random field. In this case, we explicitly 
solve the self-consistent equations for Bethe lattices with coordination numbers z = 2 and 3. 
However, this case is non-generic. For small strength of disorder A, the magnetization jumps 
from —1 to +1 at some value of the field, but for larger disorder, when the system shows 
finite avalanches, there is no jump in magnetization and the distribution function decays 



exponentially for large s. In section [V], we analyse the self-consistent equations to determine 
the form of the avalanche distribution for some other unimodal continuous distributions of 
the random field. We find that in each case for coordination number z > 4, the magne- 
tization shows a first order jump discontinuity as a function of the applied field at some 
field-strength hdisci for weak disorder. Just below h = h^sc, the avalanche distribution has 



a universal (—3/2) power-law tail. Section [VI] contains a discussion of our results, and some 
concluding remarks. Some algebraic details of the analytical solution for the rectangular 
distribution of quenched fields are relegated to two appendices. 



II. DEFINITION OF THE MODEL 

We consider a uniform Cayley tree of n generations where each non-boundary site has 
a coordination number z (see Fig. The first generation consists of a single vertex. The 
r-th generation has z(z — l) r ~ 2 vertices for r > 2. 



r = 1 

r = 2 
r = 3 
r = 4 

FIG. 2. A Cayley tree of coordination number 3 and 4 generations. 




The RFIM on this graph is defined as follows : At each vertex there is a Ising spin 
Si — ±1 which interacts with nearest neighbors through a ferromagnetic interaction J. 
There are quenched random fields hi at each site i drawn independently from a continuous 
distribution p(hi). The entire system is placed in an externally applied uniform field h. The 
Hamiltonian of the system is 
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H = -J ^2 SiSj -^h^ - hJ2 s i (!) 

<i,j> i i 

We consider the response of this system when the external field h is slowly increased 
from — oo to +00. We assume the dynamics to be zero-temperature single- spin-flip Glauber 
dynamics, i.e. a spin flip is allowed only if the process lowers energy. We assume that if 
the spin-flip is allowed, it occurs with a rate T, which is much larger than the rate at which 
the magnetic field h is increased. Thus we assume that all flippable spins relax instantly, so 
that the spin s, is always parallel to the net local field £j at the site: 

z 

Si = sign(£i) = sign(J'^2 Sj + hi + h) (2) 

i=i 

We start with h = —00, when all spins are down and slowly increase h. As we increase 
h, some sites where the quenched random field is large positive will find the net local field 
positive, and will flip up. Flipping a spin makes the local field at neighboring sites increase, 
and in turn may cause them to flip. Thus, the spins flip in clusters of variable sizes. If 
increasing h by a very small amount causes s spins to flip up together, we shall call this 
event an avalanche of size s. As the applied field increases, more and more spins flip up 
until eventually all spins are up, and further increase in h has no effect. 

III. THE SELF-CONSISTENT EQUATIONS 

The special property of the ferromagnetic RFIM that makes the analytical treatment 
possible is this: Suppose we start with h = —00, and all spins down at t — 0. Now we 
change the field slowly with time, in such a way that h(t) < h(T), for all times t < T. Then 
the configuration of spins at the final instant t = T does not depend on the detailed time 
dependence of h(t), and is the same for all histories, so long as the condition h(t) < h(T) 
for all earlier times is obeyed. In particular, if the maximum value h(T) of the field was 
reached at an earlier time t\, then the configuration at time T is exactly the same as that at 
time t\. This property is called the return point memory We may choose to increase the 
field suddenly from —00 to h(T) in a single step. Then, once the field becomes h = h(T), 
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several spins would have positive local fields. Suppose there are two or more such flippable 
sites. Then flipping any one of them up can only increase the local field at other unstable 
sites, as all couplings are ferromagnetic. Thus to reach a stable configuration, all such spins 
have to be flipped, and the final stable configuration reached is the same, and independent 
of the order in which various spins are relaxed. This property will be called the abelian 
property of relaxation. Using the symmetry between up and down spins, it is easy to see 
that the abelian property also holds whether the new value of field h" is greater or less than 
its initial value h' so long as one considers transition from a stable configuration at h' to a 
stable configuration at h" . 



We first briefly recapitulate the argument of our earlier paper [pl| which uses the abelian 
nature of spin-flips to determine the mean magnetization for any field h in the lower half of 
the hysteresis loop by setting up a self-consistent equation. 

Since the spins can be relaxed in any order, we relax them in this: first all the spins at 
generation n ( the leaf nodes) are relaxed. Then spins at generation n—1 are examined, and 
if any has a positive local field, it is flipped. Then we examine the spins at generation n — 2, 
and so on. If any spin is flipped, its descendant are reexamined for possible flips [|20|| . In 
this process, clearly the flippings of different spins of the same generation r are independent 
events. 

Suppose we pick a site at random in the tree away from the boundary, the probability 
that the local field at this site is positive, given that exactly m of its neighbors are up, is 
precisely the probability that the local field hi at this site exceeds \{z — 2m) J — h}. We 
denote this probability by p m (h). Clearly, 

POO 

Pmih) = / p(hi)dhi (3) 

J(z-2m)J-h 

Let p( r \h) be the probability that a spin on the n — r-th generation will be flipped when 
its parent spin at generation n — r — 1 is kept down, the external field is h, and each of 
its descendent spins has been relaxed. As each of the z — 1 direct descendents of a spin is 
independently up with probability p( r_1 ) ; it is straightforward to write down a recursion 
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relation for P^ in terms of P( r 1 '. For r >> 1, these probabilities tend to limiting value 



P*, which satisfies the equation |L9| 

z-l 



nh) = T, i*\p*{hT[i-p*{h)r i -™ Pm {h) (4) 



m=0 \ 111 



For the spin at O, there are z downward neighbors, and the probability that it is up is 
given by 



Prob(s = +l\h)=J2( Z ) [P\h)} m [l - P*{h)Y- m p m (h) (5) 



Because all spins deep inside the tree are equivalent, Prob(so = +1| h) determines the 
average magnetization for all sites deep inside the tree. Using Eqs. (4-5), we can determine 
the magnetization for any value of the external field h. This determines the lower half of 
the hysteresis loop. The upper half is obtained similarly. 

Now consider the state of the system at external field h, and all the flippable sites have 
been flipped. We increase the field by a small amount dh till one more site becomes unstable. 
We would like to calculate the probability that this would cause an 'avalanche' of n spin 
flips. Since all sites deep inside are equivalent, we may assume the new susceptible site is 
the site O. 

It is easy to see that this avalanche propagation is somewhat like propagation of infection 
in the contact process on the Bethe lattice. The 'infection' travels downwards from the site 
O which acts as the initiator of infection. If any site is infected, then it can cause infection 
of some of its descendents. If the descendent spin is already up, it cannot be flipped; such 
sites act as immune sites for the infection process. If the descendent spin is down, it can 
catch infection with a finite probability. Furthermore, this probability does not depend 
on whether the other 'sibling' sites catch infection. Infection of two or more descendents 
of an infected site are uncorrelated events. Thus, we can expect to find the distribution of 
avalanches on the Bethe lattice, as for the size distribution of percolation clusters on a Bethe 
lattice ETJ . However, a precise description in terms of the contact process is complicated, 



as here the infection spreads in a correlated background of 'immune' (already up ) spins, 



and the probability that a site catches infection does depend on the number of its neighbors 
that are already up. 



Y 




FIG. 3. A sub-tree Tx formed by X and its descendents. The sub-tree is rooted at X and Y is 
the parent spin of X. 



We start with the initial configuration of all spins down. Now increase the external field 
to the value h. Consider a site X at some generation r > 1 of the Cayley tree [Fig. ||. We 
call the subtree formed by X and its descendents Tx, the subtree rooted at X. We keep its 
parent spin Y at generation r — 1 down, and relax all the sites in Tx at the uniform field 
h. If X is far away from the boundary, the probability that spin at X is up is P*(h). The 
conditional probability that spin at a descendant of X is up, given that the spin at X is 
down is also P*(h). We measure the response of Tx to external perturbation by forcibly 
flipping the spin at Y ( whatever the local field there) and see how many spins in this subtree 
flip in response to this perturbation. Let Q n be the probability that the spin at X was down 
when Y was down and n spins on the subtree Tx flip up if Sy is flipped up. Here allowed 
values of n are 0, 1, 2, . . .. Clearly, we have 

oo 

P* + E Qn = 1 (6) 

71=0 

We define 
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Q{X) = £ QnX n (7) 
n=0 

Clearly, Q(x — 0) = Qo and Q(x = 1) = 1 — P*. It is straight forward to write the self- 
consistent equation for Q(x). Let us first relax all spins on T x keeping X and Y down. The 
probability that exactly m the descendents of X are turned up in this process be denoted 
by Pr{m). Clearly 

Pr(m) ={ Z ~ l \ P* m (l - py- 1 -™ (8) 
\ m J 

For a given m, the conditional probability that local field at X is such that spin remains 
down, even if Y is turned up is 1 — p m +i- Summing over m, and using the expression for 
Pr(m) above, we get 

Qo = E V^-^-^tl-Pm+l] (9) 



m=0 



m 



We can write down an expression for Qi similarly. In this case, if m of the direct 
descendents of X are up when Y is down, the local field at all the remaining z — 1 — m direct 
descendents must be such that they remain down even if X is flipped up. This probability 
is ( z ^)P* m Ql~ l ~ m - The local quenched field at X must satisfy (z - 2m) J - h > h x > 
(z — 2m — 2) J — h. The probability for this to occur is p m+ i — p m . Hence we get 

Qi = E ( Pm+ 1 - Pm) I" ~ l \ Ql~ X ~ m (10) 
m=0 V 7X1 J 

The equation determines Q n for higher n can be written down similarly. It only involves 
the probabilities Q m with m < n for the descendent spins. These recursion equations are 
expressed more simply in terms of the generating function Q(x). It is easily checked that 
the self-consistent equation for Q(x) is 

Q{X) = Q( X = 0) + X E ( Z ~ X ) (Pm + 1 - PJP*™ Q(xf- l - m (11) 



m=0 



m 



This is a polynomial equation in Q(x) of degree z — 1, whose coefficients are functions of 
h through P*(h) and p m (h). It is easily checked that for x — 1, the ansatz Q(x — 1) — 1 — P* 
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satisfies the equation, as it should. To determine Q(x) for any given external field h, we have 
to first solve the self-consistent equation for P* [ Eq. This then determines Q(x = 0) 
using Eq. [|, and then, given P* and Q(0), we solve for Q(x) by solving the (z — l)-th degree 
polynomial equation Eq. [□]. 

Finally, we express the relative frequency of avalanches of various sizes when the external 
field is increased from h to h + dh in terms of Q(x). Let G s (h)dh be the probability that 
avalanche of size s is initiated at O. We also define the generating function G(x\h) as 

oo 

G(x\h) = Y / G s (h)x s (12) 

s=l 

Consider first the calculation of G s (h) for s = 1. Let the number of descendents of O 
that are up at field h be m. For the spin at site O to be down at h , but flip up at h + dh, 
the local field ho must satisfy [(z — 2m) J — (h + dh)] < ho < [(z — 2m) J — h}. This occurs 
with probability p(z J — 2mJ — h)dh. Each of the (z — m) down neighbors of O must not flip 
up, even when so flips up. The conditional probability of this event is Qo~ m . Multiplying 
by the probability that m neighbors are up, we finally get 

G 1 (h) = J2 ( Z J p * m Qo z ~ m PizJ - 2mJ - h) (13) 

m=0 \ m J 

Arguing similarly, we can write the equation for G s (h) for s = 2, 3 etc. These equations 
simplify considerably when expressed in terms of the generating function G(x\h), and we 
get 

G(x\h) =xY J {) p * m Q(x) Z ~ m P{zJ - 2mJ - h) (14) 

m=0 \ m ) 

In numerical simulations, and experiments, it is much easier to measure the avalanche 
distribution integrated over the full hysteresis loop. To get the probability that an avalanche 
of size s will be initiated at any given site O in the interval when the external field is increased 
from hi to h 2 , we just have to integrate G(x\h) in this range. For any h, the value of dG/dx 
at x = 1 is proportional to the mean size of an avalanche, and thus to the average slope of 
the hysteresis loop at that h. 
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IV. EXPLICIT CALCULATION FOR THE RECTANGULAR DISTRIBUTION 



While the general formalism described in the previous section can be used for any distri- 
bution, and any coordination number, to calculate the avalanche distributions explicitly, we 
have to choose some specific form for the probability distribution function. In this section, 
we shall consider the specific choice of a rectangular distribution : The quenched random 
field is uniformly distributed between —A and A, so that 

p(hi) = ^ , for - A < hi < A (15) 

In this case, the cumulative probabilities p m {h) become piecewise linear functions of h, 
and /i-dependence of the distribution is easier to work out explicitly. We shall work out the 
distributions for the linear chain ( z — 2), and the 3-coordinated Bethe lattice. 



A. The Linear Chain 

The simplest illustration is for a linear chain. In this case the self-consistent equation, 
for the probability P* [ Eq. ^ ] becomes a linear equation. This is easily solved, and explicit 
expressions for Q , and Q(x) are obtained (see Appendix |S|). The different regimes showing 
different qualitative behavior of the hysteresis loops are shown in Fig. |] 

8 
4 

A 
-4 

" 8 2 4 6 8 10 

A/J 

FIG. 4. Behavior of RFIM in the magnetic field - disorder {h — A) plane for a linear chain. The 
regions A-D correspond to qualitatively different responses. In region A all spins are down and in 
region D all are up. The avalanches of finite size occur in region B and C. 
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For h < 2J — A (region A), all the spin remain down. For h > A, all spins are up (region 
D). For A < J, we get a rectangular loop and the magnetization jumps discontinuously from 
— 1 to +1 in a single infinite avalanche, and we directly go from region A to D as the field 
is increased. For A > J, we get nontrivial hysteresis loops. 

The hysteresis loops for different values of A = 0.5, 1.5 and 2.5 are shown in Fig. [5]. If 
A is sufficiently large (A > J), we find that the mean magnetization is a precisely linear 
function of the external field for a range of values of the external field h (region B in Fig. 
2). For larger h values, the magnetization shows saturation effects, and is no longer linear ( 
region C). 



M M 




(c) 

FIG. 5. Hysteresis loops for the linear chain for the rectangular distribution of quenched fields 
with different widths (a) A/ J = 0.5, (b) A/ J = 1.5 and (c) A/ J = 2.5 

The explicit forms of the generating function Q(x) are given in the Appendix |Aj. We find 
that in region B, the function Q(x) is independent of the applied field h. The distribution 
function G s (h) has a simple dependence on s of the form 

G.(h)=A 1 s(£)\ (16) 

where A\ is a constant, that depends only on J/A, and does not depend on s or h 
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A,- — ( 1 ~ J / A ) 2 (U) 
Al ~ 2A (J/A) (17) 

In region C, the mean magnetization is a nonlinear function of h. But Q(x) is still a 
rational function of x. From the explicit functional form of Q(x) and G(x\h) are given in 
the appendix [A], we find that G s (h) is of the form 

G s {h) = [A[s + A' 2 ] , for 8 > 2. (18) 

Here A[ and A' 2 have no dependence on s but are explicit functions of h. 

Integrating over h from — oo to oo we get the integrated avalanche distribution D s , 



D s = / G a (h)dh (19) 



It is easy to see from above that the integrated distribution D s also has the form 

D s = [A 2 s + B 2 ] M s>2 (20) 



.A 

where the explicit forms of the coefficients A 2 and B 2 are given in the Appendix |A|. 



B. The Case z = 3 

The analysis for the case z = 3 is very similar to the linear case. In this case, the 
self-consistent equation, for P*(h) [ Eq. | ] becomes a quadratic equation. The qualitative 
behavior of solution is very similar to the earlier case. Some details are given in Appendix ||. 
We again get regions A-D as before, but the boundaries are shifted a bit, and are shown in 
Fig. |6]. As before, in region B, the average magnetization is a linear function of h, and the 
avalanche distribution is independent of h. 
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A/ J 
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FIG. 6. Behavior of RFIM in the magnetic field - disorder (h — A) plane for Bethe lattice of 
coordination number 3. The qualitative behavior in different regions A-D is similar to that of a 
linear chain (Fig. |j). 



We find that in regime B, the distribution of avalanche sizes is given by 



G a (h) = N 



(2s)! 



>-l)!(s + 2)! 
where N is a normalization constant given by 



l-J/A)' 



It is easy to see that for large s, G s (h) varies as 



Ge ~ S 2ft S 



where 



(21) 



(22) 



(23) 



K = 4(1 -J/ A) (J/ A) 



(24) 



In region B, J/ A is always less than 1/3, and so this function always has an exponential 
decay for large s. 

In the region C, we find that the avalanche distribution is of the form 

(2a)! 



G s (h)=N' 



fs-l)!(s + 2)! 



(25) 
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where N' is a normalization constant independent of s, and Hi IS cL cL cubic polynomial in the 
external field h: 



As k is not a very simple function of h, explicit expressions for the integrated distribution 
D s are hard to write down. 



The analysis of the previous section can, in principle, be extended to higher coordination 
numbers, and other distributions of random fields. However, the self-consistent equations 
become cubic, or higher order polynomials. In principle, an explicit solution is possible for 
z < 5, but it is not very instructive. However, the qualitative behavior of solutions is easy to 
determine, and is the same for all z > 4. We shall take z = 4 in the following for simplicity. 
Since we only study the general features of the self-consistent equations, we need not pick 
a specific form for the continuous distributions of random field distribution p(hi). We shall 
only assume that it has a single maximum around zero and asymptotically go to zero at 
±oo. 

For small width ( A ) of the random field distribution i.e. for weak disorder the mag- 
netization shows a jump discontinuity as a function of the external uniform field , which 
disappears for a larger values of A fil^ . For fields h just lower than the value where the 
jump discontinuity occurs, the slope of the hysteresis curves is large, and tends to infinity as 
the field tends to the value at which the jump occurs. This indicates that large avalanches 
are more likely just before the first order jump in magnetization. 



1 



-[{9 - 53(J/A) + 119(J/A) 2 - 107(J/A) 3 } 
+ {-5 + 10(J/A) + 11(J/A) 2 } {h/A) + {3 - 9(J/A) 2 } (h/A) 2 + (h/A) 3 ] (26) 



K = 




V. GENERAL DISTRIBUTIONS 
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A/J = 2.5 : A/J = 1.5 




h/J 

FIG. 7. Magnetization as a function of increasing field for the Bethe lattice with z = 4 and the 
random field distribution given by Eq. 



For z — 4, the self-consistent equation for P*{h) [ Eq. [| ] is cubic 

aP* 3 + bP* 2 + cP* + d = (27) 

where a, b, c and c? are functions of the external field h, expressible in terms of the cumulative 
probabilities p i: i = to 3, 

a = p 3 - 3p 2 + 3p! - po 
b = 3p 2 - Qpi + 3p 
c = 3pi- 3p - 1 
d = Po 

This equation will have 1 or 3 real roots, which will vary with h. We have shown this 
variation for the real roots which lie between and 1 in Fig. || for the case where p(hi) is a 
simple distribution 

V{hi) = ^sech^ht/A) (28) 

We have also solved numerically the self-consistent equation for P* for other choices of 
p(hi), like the gaussian distribution, and for higher z(= 4,5,6). In each case we find that 
the qualitative behavior of the solution is very similar. Note that the rectangular distribution 
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discussed in the previous section is very atypical in that both the coefficients a and b vanish 
for an entire range of values of h. 

In the generic case, we find two qualitatively different behaviors: For larger values of A, 
there is only one real root for any h . For A sufficiently small, we find a range of h where 
there are 3 real solutions. There is a critical value A c of the width which separates these 
two behaviors. For the particular distribution chosen A c ~ 2.10382. 

In the first case, the real root is a continuous function of h, and correspondingly, the 
magnetization is a continuous function of h. This is the case corresponding to A = 2.5 in 
Fig. §. 

For smaller A < A c , for large ±h there is only one root , but in the intermediate region 
there are three roots. The typical variation is shown for A = 1.5 in Fig. |8|. In the increasing 
field the probability P*(h) initially takes the smallest root. As h increases , at a value 
h = hdisc , the middle and the lower roots become equal and after that both disappear 
from the real plane . At h — hdi SC the probability P*(h) jumps to the upper root . Thus 
for A < A c there is a discontinuity in P*(h) which gives rise to a first order jump in the 
magnetization curve . 



P* 




-0.5 0.5 1 1.5 2 



h/J 

FIG. 8. Variation of P*(h) with h for the Bethe lattice with z = 4, and the random field 
distribution given by Eq. ^8|. 



The field hdisc where the discontinuity of magnetization occurs, is determined by the 
condition that for this value of h, the cubic equation [ Eq. ^7| ] has two equal roots. The 
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value of P* at this point, denoted by P£ isc , satisfies the equation 

3ao^? c + 26oJ^ c + co = (29) 

where do, &o and Co are the values of a, b and c at h — hdi SC - 

We now determine the behavior of the avalanche generating function G s (h) for large s 
and h near hdisc- The behavior for large s corresponds to x near 1. So we write x = 1 — S, 
with 5 small, and h = hdisc — £■ Near hdisc, a,b, . . . vary linearly with e and 

P*w^-aVc + 0(e) (30) 

where a is a numerical constant. 

Since Q(x = 1) = 1 — P*(h) , if x differs slightly from unity Q(x) also differs from 
1 — P*(h) by a small amount. Substituting x = 1 — 5 and Q(x = 1 — 5) = 1 — P* — F(e, S) in 
the self-consistent equation for Q(x) [Eq. [H| , where both 5 and F are small, using Eq. |2^, 
we get to lowest order in 5, e and F 

F 2 + fi\f~eF — 7 2 5 = (31) 

where /3 and 7 are some constants. Thus, to lowest orders in e and 5, F is given by 



F = (l/2)[V/3 2 e + 4 7 2 5-/Vi] (32) 

Thus has leading square root singularity at x = 1 + |-f ■ Consequently, G(a;|/i) will 
also show a square root singularity £ = 1 + |-|. This implies that the Taylor expansion 
coefficients G s (h) vary as 



G s (/i) ~ s-i ( 1 + ) , for large s. (33) 



At e = 0, we get 

G s (hdisc) ~ (34) 



3 



Thus at h — hdi SC the avalanche distribution has a power law tail. 
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To calculate the integrated distribution D S} we have to integrate Eq. |33| over a range of 

2 

e values. For large s, only e < -3g- contributes significantly to the integral, and thus we get 

p s 

D s ~ ; for large s. (35) 

Thus the integrated distribution shows a robust (—5/2) power law for a range of disorder 
strength A. 



VI. DISCUSSION 

In this paper, we set up exact self-consistent equations for the avalanche distribution 
function for the RFIM on a Bethe lattice. We were able to solve these equations explicitly 
for the rectangular distribution of the quenched field, for the linear chain z = 2, and the 
3-coordinated Bethe lattice. For more general coordination numbers, and general contin- 
uous distributions of random fields, we argued that for very large disorder, the avalanche 
distribution is exponentially damped, but for small disorder, generically, one gets a jump 
in magnetization, accompanied by a square-root singularity. For field-strengths just below 
corresponding to the jump discontinuity, the avalanche distribution function has a power-law 
tail of the form s~ 3 / 2 . The integrated avalanche distribution then varies as s _5//2 for large s. 

Some unexpected features of the solution deserve mention. Firstly, we find that the 
behavior of the self-consistent equations for z = 3 is qualitatively different from that for 
z > 3. The behavior for the linear chain [z = 2) is, of course, expected to be different from 
higher z. One usually finds same behavior for all z > 2. Mathematically, the reason for this 
unusual dependence is that the mechanism of two real solutions of the polynomial equation 
merging, and both becoming unphysical (complex) is not available for z = 3. Here the self- 
consistency equation is a quadratic, and from physical arguments, at least one of the roots 
must be real. That a Bethe lattice may show non-generic behavior for low coordination 
numbers has been noted earlier by Ananikyan et al in their study of the Blume-Emery- 
Griffiths model on a Bethe lattice. These authors observed that the qualitative behavior for 



z < 6 is different from that for z > 6 22 
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The second point we want to emphasize is that here we find that the power-law tail in 
the distribution function is accompanied by the first-order jump in magnetization. Usually, 
one thinks of critical behavior and first-order transitions as mutually exclusive, as first-order 
jump pre-empts a build-up of long-ranged correlations, and all correlations remain finite- 
ranged across a first -order transition. This is clearly not the case here. In fact, the power-law 
tail in the avalanche distribution disappears, when the jump disappears. A similar situation 
occurs in equilibrium statistical mechanics in the case of a Heisenberg ferromagnet below 
the critical temperature. As the external field h is varied across zero, the magnetization 
shows a jump discontinuity, but in addition has a cusp singularity for small fields |2^| . But 
in this case the power-law tail is seen on both sides of the transition. 

Note that for most values of disorder, and the external field, the avalanche distribution 
is exponentially damped. We get robust power law tails in the distribution, only if we 
integrate the distribution over the hysteresis cycle across the magnetization jump. But, in 
this case, the control parameter h is swept across a range of values, in particular across 
a (non-equilibrium) phase transition point! In this sense, while no explicit fine-tuning is 
involved in an experimental setup, this is not a self-organized critical system in the usual 
sense of the word. Recently Pazmandi et al have argued that the hysteretic response of the 
Sherrington-Kirkpatrick model to external fields at zero temperature shows self-organized 



criticality for all values of the field [p4}| . However, this seems to be because of the presence 
of infinite-ranged interactions in that model. 

The treatment of this paper may be extended to the site-dilution case discussed by 



Tadic From the structural stability of the mechanism which leads to the cusp sin- 



gularity just before the jump-discontinuity in magnetization, it is clear that in our model, 
introduction of site dilution would not change the qualitative behavior of solutions. 

A general question concerns the behavior of the avalanches for more general probability 
distributions. Clearly, if p(hi) has a discrete part, it would give rise to jumps in pi as a 
function of h, and hence give rise to several jumps in the hysteresis loop. These could 
preempt the cusp singularity mechanism which is responsible for the power-law tails. If 
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the distribution p(hi) is continuous, but multimodal, then it is possible to have more than 



one first order jump in the magnetization ||25||. This is confirmed by explicit calculation in 



some simple cases. If p(hi) has power-law singularities, these would also lead to power-law 
singularities in pi, and hence in P*(h). Even for purely continuous distributions, the merging 
of two roots as the magnetic field varies need not always occur. For example, it is easy to 
check that for the rectangular distribution, even for z > 4, we do not get a power law tail 
for any value of A. The precise conditions necessary for the occurrence of the power-law tail 
are not yet clear to us. 

Finally, we would like to mention some open questions. Our analysis relied heavily on 
the fact that initial state was all spins down. Of course, we can start with other initial 
conditions. It would be interesting to set up self-consistent field equations for them. In 
particular, the behavior the return loop, when the external field is increased from — oo to 
some value hi, and then decreased to a lower value h,2 seems an interesting quantity to 
determine. Another extension would be to make the rate of field-sweep comparable to the 
single-spin flip rate (still assuming T=0 dynamics). This would mean some large avalanches 
in different parts of the sample could be evolving simultaneously. Then one could study 
the sweep-rate dependence of the hysteresis loops, and the frequency dependence of the 
Barkhausen noise spectra. This is perhaps of some relevance in real experimental data, 
and would also make contact with other treatments of Barkhausen noise that focus on the 
domain wall motion. 

We thank M. Barma and N. Trivedi for their useful comments on the manuscript. DD 
would like to thank the Physics Department of North Eastern Hill University, for hospitality 
during a visit there. 

APPENDIX A: 

For the case of a linear chain, the self-consistent equation, for the probability P* [ Eq. [| ] is 
a linear equation, whose solution is , 
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For h < 2 J — A, po is zero, and hence P*(h) is zero, and all the spin remain down (region 
A in Fig. |). 

For h > 2J — A, and A < J, p\ is 1 whenever po is nonzero. Then from Eq. |Al], P*(h) 
becomes 1. Thus, for A < J, we get a rectangular loop and the system changes from all 
spins down to all spins up state in a single big avalanche. 

For A > J, pi — po equals J / A and is independent of h, in the range 2J — A < h < A. Thus 
P*(h) is a linear function of h in this range, increasing from to 1. 
Defining 

e =2 (1+ A-A» (A2) 

we obtain the expression for P* as 

f for e < 



P*(/i) 



for < e < 1 - J/A (A3) 
1 for e > 1 - J/A 



Using Eq. (^), the expression for Q is, 

Qo = (1-Pi)-(P2-Pi)P*(/>) (A4) 
The generating function Q(x) obtained from the self-consistent equation [Eq.|TT|j is, 

and the generating function G(x\h) given by Eq. [TJ| becomes , 

G(x\h) = x {[Q{x)] 2 p{2J -h) + 2P*[Q(x)]p(-h) + P* 2 p(-2J - h)) (A6) 
Now if A > 2J, and -A + 2J<h<A-2J( region B in Fig. | ), 

(P2 - Pi) = (Pi - Po) = J/A, 
p(2J - mJ - /i) = — for m = 0,l,2 
and P* + Q = 1 - J/A 
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Thus 



and 



Qf x)= -J» 

^ [ ' 1 - ( J/A)x 



(A7) 



(A8) 



Expanding G(x\h) in powers of x, we get the probability distribution of avalanches in region 
B given by Eq. [L6| of sec. IV A . 

In the region C, p 2 saturates to value 1, p(—2J — h) becomes zero and (p 2 — pi) becomes 
(1 - J I A - e) . Thus we get , 

(1 - J/ A - ef 



Qo 



;i - j i a) 



In terms of P* and Qo we get 



Q(x) 



Qo + xP*[l - 2(J/A) - e] 
1 - (J/A)x 



and 



(A9) 



(A10) 



G(x\h) = ^{[P* + Q(x)] 2 -P* 2 } 



Expanding G(x\h) in powers of x we get , in region C 

Gi(h) = ^ K {(P* + Q )*-P^ 



(All) 



(A12) 



and 



G s (h) = [A[s + A' 2 



J 



for s > 2. 



2J VA, 

Here A 2 and B 2 have no dependence on s but are explicit functions of h 



A' — — — 
Al ~2A 



A' — — 
~ A) 



^(p* + Qo) 2 + ^(p- + e„)p-d-^-/VA) 



1 



2J 



2TW ( ^ + g ° )ni -?- VA) + WA) 



(A13) 



3 P^(l- 2 -l-h/Af 
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Integrating over h from — oo to oo we get the integrated avalanche distribution D Sl 



/CO 
G s {h)dh 
-oo 



(A14) 



where 



1 



(i - j/ Ay 

and, for s > 2, 



46/JA 3 47 / J_\ 4 _ 9 fJ_- 7 
Y V A J + ~6 V A ) " 5 V A 



(A15) 



(A16) 



with 



^2 
5 2 



1 



30(J/A) 
1 

15(1 - J/A) 



30 - 1 11) ( — 



135, I) 2 - 54 (k 3 



5-10.£ + 4(| 



APPENDIX B: 



For z = 3, the self-consistent equation, for P*(h) [ Eq. |] ] is a quadratic equation, 



[(P2 - Pi) - (pi - Po)]P*(h) 2 + [2(p! - po) - l]P*(h) +p = 0. 



(Bl) 



For the rectangular distribution, the coefficient of P* 2 is zero for a range of /i-values, and 
P*(h) is still a piecewise linear function of h 

' for e < 

= TWJA) for < 6 < 1 - 2(J/A) (B2) 
1 for e > 1 - 2( J/A) 

where e is defined as , 



1 h 3J N 



(B3) 



The self-consistent equation for [ Eq. |TT] ] becomes, 
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x(pi - Po) [Q{x)f + [2xP*(p 2 - p x ) - 1] Q{x) + xP* 2 (p 3 - Pa) + Qo = (B4) 

where Qo is obtained [ Eq. |9] ] as 

Qo = (1 - Pi) - 2(p 2 - Pl )P* + [(pa - Pl ) - (p 3 - p 2 )] P* 2 (B5) 

and the expression ([14] ) for G(x|/i) becomes , 

G(x\h) = x{[(Q{x)} 3 p{3J -h) + 3[Q(x)} 2 P*p(J - h) 

+3[Q(x)]P* 2 p(-J -h) + P* 3 p(-3J - h)} (B6) 

Now in the region B , 

(P3 ~Pa) = (Pa - Pi) = (Pi - Po) = 
p(3 J — 2mJ — h) = for m = to 3 
and P* + Qo = 1 - J/ A 

Solving Eq. IB3 and choosing the root which is well behaved for x near 0, we get 



1 - a/1 - 4(J/A)x(P* + Qo) 
and the expression for the integrated distribution ( IBB ) becomes 



G(x\h) = — [P* + Q(*)] 3 (B8) 



Expanding G(sc) in power series of x, we obtain the Eq. |2l] of sec. |IVB . 

In the region C, p^ saturates to the value 1, p(— 3J — h) becomes zero and (p^ — p 2 ) is 
no longer independent of h. Substituting the appropriate expressions, we find that 



and 



1 - Jl - 4(J/A)x[(l - 3(J/A) - e) + (P* + Qo)] 

Q(.) = —± 2(J /A)1 - P * 



G{x\h) = -? K {[P* + Q{x))*-F*} (BIO) 
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We note that the term inside the radical sign in Q(x), and also in G(x\h), is a simple 
linear function of x. It is thus straightforward to expand it in powers of x using binomial 
expansion. This gives us the Eq. ^5] of sec. [IV B . 
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